arXiv: 1507.05168vl [cond-mat.mtrl-sci] 18Jul2015 


Comp. Mech. manuscript No. 

(will be inserted by the editor) 


A consistent interface element 
formulation for geometrical 
and material nonlinearities 

J. Reinoso • M. Paggi 


Received: date / Accepted: date 


Abstract Decohesion undergoing large displacements 
takes place in a wide range of applications. In these 
problems, interface element formulations for large dis¬ 
placements should be used to accurately deal with cou¬ 
pled material and geometrical nonlinearities. The present 
work proposes a consistent derivation of a new inter¬ 
face element for large deformation analyses. The re¬ 
sulting compact derivation leads to a operational for¬ 
mulation that enables the accommodation of any order 
of kinematic interpolation and constitutive behavior of 
the interface. The derived interface element has been 
implemented into the finite element codes FEAP and 
ABAQUS by means of user-defined routines. The in¬ 
terplay between geometrical and material nonlinearities 
is investigated by considering two different constitutive 
models for the interface (tension cut-off and polynomial 
cohesive zone models) and small or finite deformation 
for the continuum. Numerical examples are proposed 
to assess the mesh independency of the new interface 
element and to demonstrate the robustness of the for¬ 
mulation. A comparison with experimental results for 
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peeling confirms the predictive capabilities of the for¬ 
mulation. 
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1 Introduction 

In recent years, cohesive zone models (CZMs) have been 
used in a variety of engineering applications concern¬ 
ing the formation of free surfaces due to the develop¬ 
ment of fracture processes. Relying on the seminal work 
of Barenblatt [T], CZMs have been massively incorpo¬ 
rated into computational frameworks, especially in the 
context of the nonlinear finite element (FE) method, 
as a consequence of two primary reasons: (i) the high 
versatility of the approach to accommodate different 
phenomenological fracture events, and (ii) the relative 
simplicity to numerically implement interface elements 
as user dehned subroutines into research and commer¬ 
cial FE codes. In this context, the basic ingredient that 
characterizes CZMs is the so-called nonlinear traction- 
displacement jump relationship which relates the co¬ 
hesive tractions to the relative opening and sliding dis¬ 
placements at the interface, where various contributions 
have been proposed, see BH5] for a wide review of for¬ 
mulations and a recent special issue on the topic. 

Applications cover several fields and range from quasi¬ 
static fracture in quasi-brittle solids ElIZ] with special 
attention to modelling snap-back instabilities during 
crack propagation 019 ], crack propagation in compos¬ 
ites , coupled thermo-mechanical applications |14F 

M: micromechanical and multi-scale analyses [iniiizi 
[m, fracture and contact at interfaces |20] , combination 
of friction and cohesive fracture at interfaces Eiiiia, 
and flaw-tolerance assessment in bio-inspired materi¬ 
als [BUS], among others. 

Specific contributions related to finite elements re¬ 
garded the study of the effect of the numerical inte¬ 
gration of interface elements [53], ill-conditioning situ¬ 
ations |5S], convergence issues jBIlI] and the use of a 
set of overlapping cohesive segments |5H] . 

In applications regarding thin structural elements 
subjected to large displacements, as, e.g., in biologi- 
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cal membranes, paper sheets, elastomers, viscoelastic 
materials for the encapsulation of solar cells, the com¬ 
plexity relies on the fact that during the simulation the 
deformed configuration cannot be approximated by the 
underformed one due to the occurrence of large dis¬ 
placements. Therefore, the computation of the inter¬ 
face gap (global or projected over a local reference ba¬ 
sis) according to the initial underformed geometry can 
lead to errors depending upon the specific applications 
and materials tested. Thus, large-displacement analyses 
require tracking of the surface separation, the relative 
rotations between the two sides of the interface and 
the simultaneous deformation of the two bodies sep¬ 
arated by the interface. A pioneering attempt to solve 
this problem is due to Ortiz and Pandolfi [53] , who sug¬ 
gested the adoption of a reference middle surface of the 
cohesive element in the current configuration to define 
a convenient (deformed) surface for the calculation of 
the normal and tangential directions to the interface. 
Nevertheless, their resulting formulation specified for a 
quadratic 3D interface element for matching tetrahedra 
and stemming from the differentiation of the cohesive 
tractions with respect to the normal unit vector to the 
middle surface led to a non-symmetric geometric stiff¬ 
ness matrix. In [30] . a 3D large displacement interface 
element was used based on the aforementioned formu¬ 
lation in order to simulate standard fracture mechanics 
tests in thin aluminum panels. In that case, a resid¬ 
ual with a rotation matrix updated along the deforma¬ 
tion process was considered, whereas its consistent lin¬ 
earization did not take into account the dependence of 
the standard B-operator with respect to the kinematic 
field. 

In [3I| . an alternative formulation for a 2D inter¬ 
face element in large displacements was proposed by in¬ 
troducing a non symmetric co-rotational reference sys¬ 
tem coincident with one of the two deformed sides of 
the interface. As also admitted by the authors, this co- 
rotational description leads to a very complex formula¬ 
tion of cumbersome implementation. Approximate 2D 
and 3D formulations with emphasis on the problem of 
interface fibrilation were recently proposed in [551IM] . 
In this instance, the kinematics of the interface element 
was assumed to be like that of 2D or 3D trusses under 
large displacements and rotations. More recently, an in¬ 
terface element in large displacements for fully coupled 
thermo-mechanical applications was proposed in [55] . 
The authors defined the CZM relation in a global refer¬ 
ence system, similarly to the method proposed in [55] . 
but did not consider a CZM relation that takes into 
account the contribution of different fracture modes. 
Although this could be an advantage to simplify the 
burn of the linearization of the residual, actually it re¬ 


quires the use of integrated formulations to deal with 
the transition from small to large displacement regimes 
as suggested in [55]. Moreover, to the present authors’ 
view, the geometrical contribution to the stiffness ma¬ 
trix was not clearly addressed in |33j . 

The objective of this paper is concerned with the 
development of a consistent interface element formula¬ 
tion for material and geometrical nonlinearities and the 
derivation of its corresponding finite element implemen¬ 
tation. The starting point of the consistent derivation is 
the analysis of the interface contribution to the Princi¬ 
ple of Virtual Work of the whole mechanical system, its 
virtual variation, discretization and then linearization. 
As shown in the next sections, the resulting derivation 
leads to a simple and compact operational formulation 
in which the geometric and the material contribution 
to the element stiffness matrix are clearly identified. In 
addition to this, one of the most appealing aspects of 
the model herein proposed relies on its versatility to 
accommodate any 2D and 3D finite element typologies 
along with any interface decohesion law, without any 
lack of generality. 

The article is organized as follows. In Section 2, the 
governing equations of the large displacement interface 
formulation and the corresponding finite element dis¬ 
cretization are established. The constitutive models for 
the bulk and for the interface used in this investigation 
are then briefly outlined in Section 3. In particular, a 
tension cut-off model and a polynomial CZM are con¬ 
sidered as two limit cases representative for very brittle 
or ductile interface performances. Section 4 addresses 
the main issues regarding the FE implementation in 
the context of the classical iterative Newton-Rapshon 
solution scheme. Section 5 presents a series of test prob¬ 
lems, applications to peeling and proves the robustness 
of the formulation and its ability to capture experimen¬ 
tal results related to peeling tests of very thin layers. 
Finally, the main conclusions are given in Section 6. 

2 Large displacement interface model and finite 
element formulation 

2.1 Variational framework 

The point of departure of the present formulation relies 
on the interface contribution to the expression of the 
Principle of Virtual Work of the whole system. Let us 
to assume two deformable bodies {i = 1, 2) in the 
reference configuration (denoted as Bulk-1 and Bulk-2 
in Figll]), which could have different constitutive rela¬ 
tions that characterize their mechanical performance. 
As customary, both bodies are subjected to the exter¬ 
nal body forces [i = 1,2). The boundary conditions 
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Fig. 1 A schematic definition of two bodies separated by a 
cohesive interface. 


applied on their boundaries are t® = t* on d^Q ^ and 
u* = u* on ^ {i = 1, 2). 

The bodies undergo a motion (p : x [0, t] —>■ 

where [0, t] is the time step interval, that maps the ref¬ 
erence material points (X € onto the current ma¬ 
terial points {x G ^), such that x = </)(X, t). The de¬ 
formation gradient of the transformation is defined as 
F dxp{^,t), with the Jacobian J = det[F] and dx 
denoting the partial derivative with respect to the refer¬ 
ence frame. Moreover, it is supposed that the interface 
between both solids is characterized by the presence of 
a cohesive surface Sq. 

Focusing our attention on the analysis of the inter¬ 
face between the solids, the contribution of the interface 
cohesive tractions T, acting on ^o, to the Principle of 
Virtual Work of the mechanical system in the reference 
configuration is: 


-^intf(gloc) — / Sioc'FdS' (1) 

JSo 


where gioc is the gap vector accounting for opening and 
sliding displacements between the two sides of the in¬ 
terface. Note that, due to the geometrical nonlinearity, 
the traction vector previously defined corresponds to 
the nominal first Piola-Kirchhoff tractions related to 
the local basis of the interface in the reference configu¬ 
ration. 

It is worth noting that in the large deformation set¬ 
ting, the gaps vector vanishes when the body undergoes 
rigid body motions, thus confirming the frame indiffer¬ 
ence of the formulation proposed in this paper. 


Reference Configuration Current Configuration 



Fig. 2 Kinematics definition of the interface along the de¬ 
formation process. 

The virtual variation of TTintf according to the prin¬ 
ciple of virtual displacements reads: 

TdS 

Js. ^ , 2 ) 

In case of large displacements, the updated coordi¬ 
nates of a generic point are given by x = X -|- u, see 
FiglSl 

As is generally proposed for interface formulations, 
it is convenient to define a middle line (in the 2D case) 
in the updated configuration by averaging the position 
vectors and the displacement fields of the upper and 
lower sides of the interface, see FiglS] after performing 
a standard discretization process. Hence, the position 
vector X of a generic point along this middle line can be 
determined by pre-multiplying the positioning vector x 
by an averaging operator M: 

X = Mx (3) 

2.2 Finite element formulation 

Based on isoparametric interpolation, the position vec¬ 
tor at the interface can be approximated through: 

X x'^ = Nx" (4) 

where x” denotes the nodal position vector (the su¬ 
perscript n identifies nodal quantities), and N is the 
the operator that collects the shape functions and it 
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Fig. 3 Sketch of the interface element with node numbering 
and integration points. 


depends on the natural coordinate of the element In¬ 
troducing now the discretization of the interface into 
Eq.Q, the interpolated average position vector yields: 

X = NMx" (5) 

Similarly, the coordinates of the points belonging to 
the middle line in the reference configuration, X, and 
their displacement vector, u, can be computed via a 
standard interpolation procedure from the nodal quan¬ 
tities 


X ^ X"" = NMX”, u ^ u® = NMd (6) 

where X” and d denote the position vector of the nodes 
in the reference configuration and their nodal displace¬ 
ment vector, respectively. 

In 2D, the tangential and the normal vectors t and 
n to the middle line of the interface element used to 
define the local frame are given by: 

dx'^ 

t = n.t = 0 (7) 

Note that in 3D application, in line with the derivation 
proposed in Ortiz and Pandolfi [52] , the convective tan¬ 
gential and normal vectors to the middle surface of the 
interface element (ti, t 2 and n) used to define the local 
frame are determined via differentiation of the average 
coordinates with respect to the natural coordinates ^ 
and 77 : 

i9x® X® 

*1 = ^’ <^2 = ^, n = tixt2 (8) 

The gap vector in the reference cartesian frame, g, 
can be obtained by pre-multiplying the nodal displace¬ 
ment vector d by a suitable operator L which pro¬ 
vides the difference between the displacements of the 


upper and the lower bodies at the interface. Accord¬ 
ingly, within the FE discretization we have: 

g - g^ = NLd (9) 

The constitutive relation for the interface, i.e., the so- 
called cohesive zone model (CZM), is usually provided 
in a local frame defined by the normal and the tangen¬ 
tial vectors to the average line of the interface element 
in order to distinguish between fracture Modes I and 
II, as introduced in Eq. ©. Therefore, the gap vector in 
this local frame, gioc, has to be computed by multiply¬ 
ing the gap vector in the reference frame by a rotation 
operator R: 

gioc = R(u)g (10) 

It is remarkable to note that, in case of large displace¬ 
ments, the operator R(u) is a function of the displace¬ 
ment field. Its expression is detailed in Section 4 for the 
2D case that represents the main scope of the present 
work (the 3D version can be derived adapting the for¬ 
mulation here developed). Consequently, a consistent 
formulation must take into account this dependency in 
the subsequent linearization of the discretized version 
of the interface contribution to the Principle of Vir¬ 
tual Work within the classical Newton-Raphson itera¬ 
tive solution scheme. This dependency will lead to the 
so-called geometric contribution to the element stiffness 
matrix. Introducing the FE discretization, Eq. (nni) can 
be rewritten as: 


gf„, = R(d)NLd 


( 11 ) 


Examining the terms entering the virtual variation 
of the virtual work in Eq.Q, the partial derivative 
(9gioc/9u) is approximated by: 


^ R(d)NL + ^^^NLd (12) 

au ad ad 

where the differentiation of the second order tensor R 
with respect to the components of the vector d leads to 
a third order tensor. Note that in (115]) the formulation 
is simplified by omitting the second derivative of the ro¬ 
tation matrix with respect to the displacement vector. 
This vanishes in case of linear displacement interpola¬ 
tion under the assumption that the norm of the tangent 
vector t does not depend on the displacement field. In 
this regard, we also assessed the role of this term based 
on representative numerical tests adopting alternative 
interpolation schemes. It was found that this term has 
an almost negligible effect on the results and it is there¬ 
fore reasonable to be neglected. 

The operator B = NL is now introduced and Eq. (HH) 
can be rephrased as: 


dg:. 


dR, 


^ = RB + ^Bd 


(13) 
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The matrices R and B are evaluated at the element 
level, though the typical superscript (e) has been omit¬ 
ted here to simplify notation. 

Inserting this intermediate result into the discretized 
version of Eq. where u is simply replaced by d, the 
following general formulation valid for any kind of in¬ 
terface element topology dealing with geometric and 
material nonlinearities is derived: 


After some algebra we obtain the following result: 


K® = / B^R^CRB dS” 

JSo 

,5RT 


-b 


/So L 


2 B' 


5d 


aRT AR_ , 

T-l-d^B^^:r^C—Bd (19) 


5d dd 


rp rp c/R rp rp C/R^ 

TB^rTc—B d + d^B^^—CRB 
ad ad 


d5 


6W 


intf 



RB 


5R 

dd 


Bd 


T 

TdS 


(14) 


Summarizing, the tangent stiffness matrix which ac¬ 
counts for both the material and the geometric contri¬ 
butions reads: 


The solution of the variational equation = 

= 0 V5d results in the equations set = 0, 
where is a nonlinear function of the unknown d and 
assumes the role of the residual vector in the Newton- 
Raphson iterative scheme: 

= ^ (RB + ^Bd) TdS (15) 


_ T^e I 

mat ' geom 


TC^ = 


L 


B^R^CRBdS' 


K" 


I So 


2 B^-r^T -b d^B^^zp^C—Bd 


dd 


rp rp dR rp rp dR^ 

B^R^C—Bd + d^B^^—CRB 

dd dd 


dd dd 
d^ 


(20a) 

(20b) 


(20c) 


which leads to the following equations set for the com¬ 
putation of the corrector Z\d at each iteration k: 


= (16a) 

dfc-si = d'= -b Ad (16b) 


To alleviate the notation, the superscript k will be omit¬ 
ted in the sequel. Following standard arguments of non¬ 
linear FE formulations, the element stiffness matrix K® 
is obtained from the linearization of the residual, i.e., 
K® = dfintf/dd and it is evaluated by using the dis¬ 
placement held solution at the iteration k: 


K" 



dR"^ 

^d“ 


T 


B^RT -bd^B^ 



dT' 

dd 


dS 


(17) 


In case of small displacements, Eq. o reduces to 
the standard form of the material contribution to the 
element stiffness matrix, Eq. (I20bl) . In case of large dis¬ 
placements, the complete tangent stiffness matrix is 
composed of four terms, see Ea. (l20cl) . Only the hrst 
one, involving the computation of the cohesive traction 
vector T is not symmetric. Therefore, a nonsymmetric 
solver has to be used. However, in case of a symmetric 
constitutive matrix C for the interface, as it happens 
in case of the same CZM parameters for Mode I and 
Mode II, we explored the possibility to neglect the non 
symmetric contribution to Kg. The examples discussed 
in Section 5 will show that this will not affect the ac¬ 
curacy of the solution and slightly penalize the con¬ 
vergence rate. The omission of this term, on the other 
hand, makes it possible the use of symmetric solvers. 
This fact represents an obvious advantage in case of 
massive computations. 


In this derivation, as it was already stated, the second- 
order differentiation of the rotation matrix R was omit¬ 
ted. 

The derivative of the cohesive traction vector T with 
respect to the displacement vector d can be determined 
via a chain rule differentiation: 


dT 


dT gioc 
dgioc dd 


RB 


dd 


Bd 


(18) 


where C = —- represents the tangent interface con- 

^Sloc 

stitutive matrix whose expression will be detailed in the 
next section. 


3 Material models 

With reference to the continuum, to assess the effect 
of large displacements on debonding of thin structures, 
both a small deformation and a large deformation ver¬ 
sions of a standard homogeneous isotropic hyperelastic 
material model are considered in the sequel. The deriva¬ 
tion is here omitted for the sake of brevity. The readers 
can refer to [36] for more details. 

Regarding the cohesive traction vector T, with the 
aim of quantifying the role of the geometric nonlinear 
effects along the decohesion process, two different types 
of interface constitutive laws are examined. 
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First, a tension cut-off CZM is considered, with un¬ 
coupled Mode I and Mode II deformation. This type 
of CZM has the advantage of allowing for closed form 
solutions for specific testing configurations, like for the 
double cantilever beam test [37]. The stiffness of the 
CZM can be related to the Young’s modulus E and to 
the thickness of the adhesive h, i.e. k = CTmax/^nc ^ E/h 
where (Tmax and l^c denote the critical traction for dam¬ 
age initiation and the critical relative displacement, re¬ 
spectively. In this approach, when the crack sliding or 
opening displacements overcome a critical value corre¬ 
sponding to the achievement of the adhesive strength, 
Cmax, the interface suddenly debonds. Since such crit¬ 
ical relative displacements for failure are very small 
quantities in applications, the process zone size is ex¬ 
pected to be quite small and limited within a region 
very close to the real crack tip, where displacements 
are moderately small. 

The cohesive traction vector T = (r, cr)'^ reads: 


—”^m.ax 


^ —^max 


5'loc,t 

^tc 

5loc,n 

^nc 


( 21 a) 

( 21 b) 


where (/ioc,t and 5 ioc,n are the tangential and normal 
components of the gap vector gioc, whereas Itc and 
/nc are the critical sliding and opening displacements. 
The tangent constitutive matrix stemming from the lin¬ 
earization of the CZM tractions with respect to the gap 
vector is: 


C = 


"Rnax 

^tc 

0 


0 


<^max 




nc 


( 22 ) 


In this case, in line with the previous arguments, the re¬ 
sulting interface element stiffness matrix is always sym¬ 
metric. 

As second formulation, we consider the polynomial 
CZM by Tvergaard [35| as an example of an interface 
constitutive relation where the cohesive traction vec¬ 
tor T = (r, cr)"^ is a nonlinear function of the sliding 
and opening displacements with a softening branch af¬ 
ter reaching the maximum cohesive tractions. For the 
same values of the parameters T^a.x, CTmax and of the 
initial stiffness as for the tension cut-off model, this 
CZM has a larger fracture energy and therefore a more 
widespread process zone is expected. In this model, the 
cohesive tractions are given by 

T =rmax^^^F(A) (23a) 

^tc 

^ =a„,ax^F(A) (23b) 


where 


P(A) 


— (1-2A + A2),for0< A< 1 
0 , otherwise 



(24a) 

(24b) 


For this CZM, the tangent constitutive matrix reads: 


C = 


P 

^tc 


‘ max . Rnax 


9\oc,t 


7 

He 

9\oc,r, 




P 


gioc.t dP dX 
^tc dX Ogioc t 

dP dX 

dX Ogioc^n 

dP dX 

dX dg\oc,t 
gioc,n dP dX 

^nc dX Ogioc,n 


(25) 


4 Matrix operators for finite element 
implementation 

This section covers the main features concerning the 
numerical implementation of the large displacement in¬ 
terface element formulation proposed in the previous 
section. According to the derivation presented in Sec¬ 
tion 2 , we restrict our attention to the implementation 
of the element in a 2D version, although Eqs.(18) and 
(19) are the same for 3D problems, provided that a mid¬ 
dle surface is introduced in analogy with the middle line 
for the 2D case. 

Let us to consider a 4 node bilinear interface ele¬ 
ment, see Fig 131 The corresponding shape functions to 

accomplish the numerical integration are A^i = —(1 — 0 

and N 2 = ^(1 + 0- Each node has two degrees of free¬ 
dom, so that the nodal position and displacement vec¬ 
tors are arranged as: 


X = (Yi, Yi, X 2 , Y 2 , X 3 , Fa, Y 4 , Y 4 )^ (26a) 

d = (ui,Ui,U2,U2,U3,U3,U4,U4)'’" (26b) 

where Y^, Yi identifies the cartesian coordinates corre¬ 
sponding to the node i and ut and Vi stands for the 
corresponding displacements along the Y and Y direc¬ 
tions. 

The gap and the traction vectors that characterizes 
the CZM are evaluated in correspondence of each inte¬ 
gration point, so that: 

Sloe — (gioC,t; gioC,n) (27a) 

T=(T,a)T (27b) 
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Next, the matrix operators defined in Section 2 to 
determine the coordinates and the gaps of the points 
belonging to the interface middle line take the form: 


N = [A^il iVsl] 

-10 0 1 
0 -I I 0 


M=- 

2 


10 0 1 
0 I I 0 


L = 


(28a) 

(28b) 


where 0 is a 2 x 2 null matrix and I is a 2 x 2 identity 
matrix. As previously stated, the CMZ is evaluated at 
the local reference system that is defined by the tan¬ 
gential vector t and the normal vector n to this middle 
line. Thus, the rotation operator yields: 


R = 


tx ty 

Tlx 


(29) 


where: 


tx - Tiy - 


ty - Tlx - 


X2 U2 U3 — — Ui — -^4 — 'U4 

m\ 

Y2 + V2 + Y3 + V3 - Yi - Vi - Y4 - V4 

2||t|| 


and the symbol ||.|| denotes the Euclidean norm of the 
corresponding vector. Note that, differing from previous 
interface formulations [52] j this operator is evaluated at 
each integration point at the element level. 

The operator stemming from the third order tensor 
has to be computed with care. For the present particu¬ 
lar case it renders: 


Data: Given: X®, d, Ad, At 
if nlgeom > 0 then 
I Update the geometry: x® = X® -|- d; 

else 

I x®=X®; 

end 

Construct L; 

Loop over the integration points; 

for j = 1 to 2 do 

Evaluate shape functions and derivatives; 

Compute the local basis vectors [t,n] and rotation 
matrix R; 

Construct N; 

Perform B = NL; 

Compute RB (RB)^; 

Evaluate T, C according to the selected CZM; 

end 

if nlgeom > 0 then 

Perform: -^Bd; 

ad ’ 

Construct the geometrical stiffness matrix 
Ea. (l20cll : 

Compute the geometric part of the residual vector 
^intf > second term of Eq.(15); 
else 

Compute the material part of the residual vector 
first term of Eq.(15); 

end 

Evaluate the material contribution to the stiffness 
matrix Ea. (l2nbll : 

if nlgeom > 0 then 

I K® = -b K®g„^ 

else 

I K® = KS,,t 

end 

Update stiffness matrix and r.h.s. vector; 

Algorithm 1: Numerical implementation of the 
large displacement interface element. 


, —a -b +a +b +a +b —a -b 

vrrBd = , , , , 30) 

i9d [—6 +a +b —a +b —a —b -|-aj ' ^ 

, —Niui — N2U2 + N2U3 + N1U4 , , 

where a = -- and b = 

2||t|| 

—Nivi - N2V2 + N2V3 + N1V4 

2p ^ • 

The remaining operations to accomplish are straight¬ 
forward, and therefore are omitted here for the sake of 
conciseness. 

The algorithmic treatment of the proposed formu¬ 
lation, implemented by the present authors both in the 
Finite Element Analysis Program FEAP and in the 
commercial software Abaqus as user defined subrou¬ 
tines, is summarized in the following sequence of main 
operations, see Algorithm 1. Note that the external loop 
over j refers to the numerical integration, whereas the 
variable nlgeom indicates a flag defined in the input file 
to select between: (i) small displacement formulation 
{nlgeom = 0) or (ii) large displacement formulation 
{nlgeom = 1). 


5 Numerical examples 

In this section, the numerical performance of the pro¬ 
posed element is illustrated. To this aim, two appli¬ 
cations are selected. First, we investigate the element 
capabilities through a series of benchmark problems to 
highlight the principal capabilities that the present ele¬ 
ment formulation incorporates. Second, a structural ap¬ 
plication consisting of a peeling test is addressed. Small 
or finite displacement formulations for the continuum 
and for the interface element are examined, along with 
two different CZM formulations. 

5.1 Benchmark tests 

A preliminary test problem shown in Fig |4(a)| is an¬ 
alyzed in order to assess the performance of the new 
interface element for large displacements as compared 
to a standard formulation for small displacement. This 
benchmark problem aims at investigating the interplay 
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between geometric and material nonlinearities, an issue 
not yet rigorously quantified in the related literature. 
Therefore, as was previously stated, we consider small 
or large deformation hyperelastic material models for 
the continuum in order to investigate the role played 
by the different interface element formulations in these 
two cases. Hence, each numerical test will be identified 
by labels X — Y{Z). Whereas X refers to the interface 
element formulation and it can be S' or L depending on 
the small or the large displacement formulation used, 
the second label Y stands for the constitutive model of 
the continuum and again it can be S or L depending 
on the small or large deformation theory adopted for 
the hyperelastic material. The symbol Z in parenthesis 
denotes the solver used (s for symmetric solver or u for 
a non symmetric solver). 

In this benchmark problem, two blocks of lateral size 

I mm and different heights are discretized by a single 
finite element each. The lower block of unit height has 
a Young’s modulus Ei = 500 MPa in order to simulate 
an almost rigid substrate, whereas the upper block has 
a height of 0.1 mm and a Young modulus E 2 = 5 MPa 
to simulate a highly deformable elastomeric tape or a 
paper sheet. Both materials have a vanishing Poisson’s 
ratio and the simulations are conducted under plane 
strain assumption. An interface element is placed be¬ 
tween the two blocks and its CZM can have either a 
tension cut-off or a polynomial form, see Fig |4(b)| for the 
Mode I relations (Section 3). The values of the CZM pa¬ 
rameters Umax and Tmax are selected the same for both 
the tension cut-off and the polynomial CZMs, to per¬ 
form a consistent comparison. On the other hand, the 
fracture energy of the tension cut-off model is much 
smaller than that of the polynomial CZM, as it can 
be readily visualized in Fig |4(b)] from the different area 
under the respective traction-separation relations. The 
value of Inc in the tension cut-off model has been se¬ 
lected so that the tension cut-off curve is tangent to 
the polynomial CZM at 510 c,n = 0. 

As far as the boundary conditions are concerned, 
a non-uniform Mixed Mode debonding problem is sim¬ 
ulated by restraining the first block at the basis and 
imposing a linear vertical displacement variation to the 
upper side of the second block. This testing configura¬ 
tion has been chosen because leads to large displace¬ 
ments at the interface during the deformation process 
and therefore Mode Mixity m- In other types of load¬ 
ing, such as uniform Mode I debonding, uniform Mode 

II debonding or uniform Mixed Mode debonding, the 
response of the interface element would be in fact the 
same regardless of the small or large displacement for¬ 
mulation used, since the orientation of the local frame 
does not change during the deformation process. These 



(a) Sketch of the test problem 



Fig. 4 Sketch of the geometry of the 2-blocks Mixed Mode 
test problem and illustration of the normal cohesive traction- 
normal gap CZM relations used in the simulations. 


trends are consistent with the results reported in [32], 
although they were based on an approximated large 
displacement formulation for the interface element de¬ 
duced by the kinematics of a beam element in large 
displacements and rotations. 


The vertical reaction force E in the constrained node 
is plotted vs. the imposed displacement A in Figj^for 
an interface tougher in Mode I than in Mode II, i.e., 
with CTmax/Tmax = 10 and Inc = he = O.I mm. The re¬ 


sults for the tension cut-off CZM are shown in Fig 5(a) 


and those for the polynomial CZM are depicted in Fig 5(b) 


Examining the curves related to the tension cut-off CZM, 
Fig 5(a) we note that, in case of the small deformation 
theory for the continuum, curves labeled (S' — S) and 
{L — S) are almost coincident before the peak load, i.e., 
before complete debonding of the first Gauss point of 
the interface. Note that the labels (s) and {u) make 
reference to the symmetric or unsymmetric character 
of the formulation. In this sense, the small-displacement 
case always leads to a symmetric stiffness matrix, whereas 
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the large displacement case provides a general unsym- 
metric formulation, in which the role of the term that 
break the symmetry (see Eq.(20c)) is examined. 

After the peak load, an abrupt reduction in the in¬ 
terface load-carrying capacity takes place and the struc¬ 
tural response is characterized by a lower stiffness until 
complete decohesion takes place. In this respect, the 
large displacement formulation for the interface pre¬ 
dicts a much lower displacement jump for complete 
debonding. A similar difference can be noticed in case of 
a large deformation formulation for the continuum, see 
curves labeled (S — L) and (L — L). Since this CZM does 
not consider coupling between Mode I and Mode II, i.e., 
the matrix C is diagonal, it makes sense to compare 
the results for the large displacements interface element 
formulation by considering the complete expression of 
the tangent stiffness matrix and using a non symmetric 
solver or the approximate symmetric expression and us¬ 
ing a symmetric solver. As it can be seen from Fig |5(a){ 
the results are coincident. 

In case of the polynomial CZM, the results have the 
same trend as for the tension cut-off, just with much 
smoother curves during the debonding process. In light 
of the previous arguments, it is worth noting that for a 
given assumption regarding the kinematics of the con¬ 
tinuum, the use of a large displacement formulation for 
the interface instead of its small displacement counter¬ 
part has a predominant effect on the softening branch 
of the F — A response. In this case, since the matrix C 
is not symmetric due to the expression of the CZM, a 
non symmetric solver has been always used and the full 
expression for the geometric stiffness matrix has been 
retained in the computations. 

Examining the element performance in case of an 
interface with the same fracture parameters in Mode 

I and in Mode II (crmax/Tmax = 1), see EiglHl we find 
that the discrepancy between the predictions in case of 
large or small interface element formulations are min¬ 
imal. On the other hand, small or large displacement 
formulations for the continuum significantly affects the 
post-peak branch. Since in this case the matrix C is 
symmetric for both the CZM formulations, the use of 
the complete non symmetric tangent stiffness matrix or 
its symmetric version by neglecting the non symmet¬ 
ric contribution to its geometric component have been 
compared and the results are again coincident. 

Finally, the last scenario to be inspected is repre¬ 
sented by the case of an interface much tougher in Mode 

II than in Mode I (o'max/'rmax = 0.1), see FiglTl As far 
as the choice of the symmetric or the non symmetric 
solver is concerned, the same comments to the case 
when Mode I prevails over Mode II apply. In this in¬ 
stance, since the loading test is predominantly in Mode 




Fig. 5 Total force vs. imposed displacement for the test 
problem in Fig |4(a)] and different CZM formulations for an in¬ 
terface tougher in Mode I than in Mode II (crmax/Tmax = 10). 
(S — S)\ small displacement interface element & small de¬ 
formation continuum; {L — S)\ large displacement interface 
element & small deformation continuum; (S — L): small dis¬ 
placement interface element & large deformation continuum; 
{L — L)\ large displacement interface element & large defor¬ 
mation continuum. 

I, we observe a much lower peeling force F in this bench¬ 
mark test. Additionally, slight discrepancies between 
the numerical predictions using large or small interface 
element formulations are noticed. 

Therefore, it is possible to draw the practical con¬ 
clusion that the large displacement formulation for the 
interface should be primarily used in case of applica¬ 
tions with CTmax > fmax, as, e.g., in fibrilation problems 
where the shear strength of cellulose or polymeric fibrils 
is almost negligible as compared to their axial strength. 

5.2 Structural application: peeling test and 
comparison with experiments 

Examining now a structural problem where the large 
displacement formulation for the interface element is 
deemed to be crucial, a peeling test where a thin layer 
is pulled from al almost rigid substrate by the action of 
a vertical displacement imposed to the top right corner 
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(b) Polynomial CZM 


Fig. 6 Total force vs. imposed displacement for the test 
problem in Fig |4(a)] and different CZM formulations for an 
interface with the same toughness in Mode 1 and in Mode 
II (u'max/'rmax = !)■ {S — S)'. Small displacement interface 
element & small deformation continuum; (L — S): large dis¬ 
placement interface element & small deformation continuum; 
{S — L)\ small displacement interface element & large de¬ 
formation continuum; {L — L)\ large displacement interface 
element & large deformation continuum. 


Fig. 7 Total force vs. imposed displacement for the test 
problem in Fig |4(a)] and different CZM formulations for an 
interface tougher in Mode If than in Mode 1 (tTmax/Tmax = 
0.1). (S — S): small displacement interface element & small 
deformation continuum; (L — S): large displacement interface 
element & small deformation continuum; (S — L): small dis¬ 
placement interface element & large deformation continuum; 
{L — L)\ large displacement interface element & large defor¬ 
mation continuum. 


is considered (see the final deformed shapes in case of 
{S — S) or {L — L) formulations in FiglS]). The mate¬ 
rial parameters for the bulks and for the CZMs (tension 
cut-off and polynomial CZMs) are the same as in the 
previous example, considering the case of an interface 
tougher in Mode I than in Mode II (cTmax/fmax = 10, 
where the large displacement formulation for the inter¬ 
face element was found to significantly differ from the 
small displacement one. A non symmetric solver and 
the complete expression for the tangent stiffness ma¬ 
trix are used. 

The force-displacement curves for different kinemat¬ 
ics formulations are compared in FigEl As a general 
trend, the large displacement formulation for the inter¬ 
face element leads to lower peak loads as compared to 
its small displacement counterpart, for a given kinemat- 
ical model of the continuum. Large differences among 
the predictions of the formulations can also be observed 
as far as the softening branches are concerned. 



DispI, y 
1.5485 


0.4 


iO.2 

to 

-0.00435 


Fig. 8 Deformed meshes of the peeling test in case of small 
displacement formulations for the continuum and the inter¬ 
face element (left, (S — S)) or large displacement formulations 
(right, (L - L)). 
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Fig. 9 Total force vs. imposed displacement for the peeling 
test in Fig[7] and different CZM formulations for an inter¬ 
face tougher in Mode I than in Mode II (crmax/Tmax = 10)). 
{S — S): small displacement interface element & small de¬ 
formation continuum; {L — S): large displacement interface 
element & small deformation continuum; (S — L): small dis¬ 
placement interface element & large deformation continuum; 
{L — L): large displacement interface element & large defor¬ 
mation continuum. 

Additionally, some comments on the convergence 
of the formulation herein proposed have to be added. 
First, from the numerical point of view, the interface 
element formulation was found to be quite stable, with 
the appearance of small oscillations in the softening 
branches for the peeling test only in case of the tension 
cut-off CZM. These effects are caused by the sharp dis¬ 
continuity in the traction-gap constitutive relation lead¬ 
ing to small jumps in the load when debonding takes 
place in a given Gauss point, see Figs |9(a)} These small 
oscillations disappear in case of the polynomial CZM, 
since a softening is included in the interface constitutive 
relation. 

Second, for the sake of completeness, the mesh con¬ 
vergence of the method is tested by considering the 
polynomial CZM and performing peeling tests as in 
Fig |9(b)| for the (L-L) case, with different mesh refine¬ 
ment for the bulk and the interface in the horizontal di¬ 
rection. In the coarsest discretization (mesh 1), only 25 



Fig. 10 Mesh convergence study for the peeling test in FiglH] 
Mesh 1 corresponds to 25 elements along the interface, mesh 
2 to 50 elements, mesh 3 to 100 elements and mesh 4 to 200 
elements. 

elements along the interface are used. Finer meshes with 
50 elements (mesh 2), 100 elements (mesh 3) and 200 el¬ 
ements (mesh 4) along the interface are also considered. 
The corresponding load-displacement curves are shown 
in Fig llOl where an excellent mesh-independency even 
for the coarsest mesh can be observed. 

The predictive capabilities of the proposed formula¬ 
tion are finally checked against experimental results. 
To this aim, a 90° peeling of a backsheet (0.1 mm 
thick, Young’s modulus 2.8 GPa, vanishing Poisson’s 
ratio, hyperelastic material) from a glass substrate (4 
mm thick, Young’s modulus 73 GPa, vanishing Pois¬ 
son’s ratio, linear elastic material) is simulated by mod¬ 
elling the adhesive response of the Epoxy Vynil Acetate 
(EVA) interlayer via the polynomial CZM used in the 
previous examples. The parameters to be identified are 
the peak cohesive traction (Tmax and the fracture en¬ 
ergy Gic, which is proportional to the critical open¬ 
ing displacement Inc- The same parameters for Mode I 
and Mode II deformation are used and the symmetric 
formulation of the interface element for large displace¬ 
ment analyses is adopted. The test is conducted under 
plane strain conditions. Mesh and boundary conditions 
are analogous to those displayed in Eig.8, although the 
upper layer is much thinner in the present problem. 
Eour finite elements are used to discretize the back- 
sheet through its thickness and 200 finite elements are 
used along the interface, considering an initial bonded 
length of 50 mm. 

For this test, essential to ascertain the reliability of 
backsheet bonding in photovoltaic systems, experimen¬ 
tal results are reported in m- Since the material pa¬ 
rameters and the exact dimensions corresponding to the 
experimental force-peel extension curve are not listed 
in m. we use values conforming to the standard ma¬ 
terials used in PV production, see [32]. Another source 
of uncertainty regards the way the peeling extension 
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Fig. 11 Peeling of a backsheet from a glass substrate: nu¬ 
merical vs. experimental results taken from un¬ 
is measured, since no details are provided in [JT]. In 
the numerical simulation we predict the peeling exten¬ 
sion via the location of the fictitious crack tip position. 
Keeping in mind that, alternatively, the position of the 
real crack tip could be used, this choice makes a cer¬ 
tain difference especially in the pre-peak branch of the 
force-peel extension curve. 

Numerical results are shown in FigfTTlfor Gjc = 5.4 
N/mm, which corresponds to the steady-state peeling 
force measured in experiments. A series of curves ob¬ 
tained by varying (Tmax in the range from 3.6 to 28.8 
N/mm are displayed and are used to identify the value 
of CTmax which provides the best agreement with experi¬ 
ments. In the present case, we found Umax 5.8 N/mm, 
corresponding to the black dashed curve in FigfTTl The 
numerical methods is able to predict the stead-state 
peeling force very well, in excellent agreement with ex¬ 
periments. The pre-peak response, which has the high¬ 
est degree of inaccuracy in the real tests, is in any case 
reasonably well reproduced. 


6 Conclusions and outlook 

In this paper, a consistent derivation of an interface 
element for large displacements applications has been 
proposed. 

The present theory finds its variational basis in the 
interface contribution to the Principle of Virtual Work 
of the whole mechanical system. Our present work, dif¬ 
fering from alternative formulations presented in the 
literature, furnishes a consistent derivation of the inter¬ 
face model involving large deformations. Particularly, 
the cohesive model herein developed takes into account 
the full finite kinematics in which the material and the 
geometrical contributions to the element stiffness ma¬ 
trix are clearly determined. 


The corresponding finite element discretization of 
the interface model has been accomplished based on a 
linear two-dimensional zero-thickness interface element 
for which the fundamental operators and the implemen¬ 
tation details have been addressed. The compact and 
consistent theoretical derivation allows its straightfor¬ 
ward generalization to different orders of the kinematic 
interpolation and to 3D topologies. 

Numerical applications using two different cohesive 
interface laws and with small or finite deformation kine¬ 
matic assumptions for the continuum have been ex¬ 
amined in order to assess the interface element per¬ 
formance. Particularly, concerning the interface consti¬ 
tutive models, two laws have been selected: (1) the so- 
called tension cut-off model, that assimilates a quasi- 
brittle behavior of the interface, and (2) the polyno¬ 
mial based Tvergaard model that was adopted for sim¬ 
ulating a ductile interface. The numerical results have 
proven the applicability of the interface element pro¬ 
posed especially in terms of its satisfactorily numerical 
convergence in achieving equilibrium solutions, along 
with a minimal mesh sensitivity. The predictive char¬ 
acter of the method has been demonstrated through 
the simulation of a peeling test of a backsheet from a 
glass substrate, in which the ability of the formulation 
to capture the nonlinear character of the experimental 
trend is noteworthy. 

In closing, we would like to emphasize that the de¬ 
veloped model is particularly promising in addressing 
real situations undergoing large displacements, which 
commonly take place in a wide range of engineering 
and biomechanical applications. This fact has been evi¬ 
denced in the peeling simulations included in this inves¬ 
tigation where the role of finite displacements has been 
highlighted. In case of problems involving thin layers 
with an interface much tougher in Mode I than in Mode 
II, as in fibrilation problems, the proposed interface el¬ 
ement for large displacements is recommended to be 
used instead of its small displacement counterpart, to 
avoid a significant overestimation of the peeling force, 
as shown in the examples discussed in the present study. 
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